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■OUTKV 


This  user's  asnusl  Is  designed  to  ssslst  the  ■sthematlclsn  or  progrsaaer 
using  the  BRL  Polygsnaa  Function  subroutine  .  FORTRAN  symbols  for  variables 
and  arithmetic  operations  are  used  In  the  body  of  the  report  for  consistency 
with  excerpts  from  the  coding  . 

As  an  aid  to  the  reader  unfamiliar  with  standard  FORTRAN,  the  following 
8]nibols  are  defined: 


Algebraic 

F(»TRAN 

Symbol 

Operation 

Notation 

Notation 

1  .  + 

add 

a  +  b 

-  A  +  B 

2  . 

subtract 

a  -  b 

-  A  -  B 

3  .  * 

multiply 

a  X  b 

-  A  *  B 

4.  / 

divide 

a  t  b 

-  A  /  B 

Numbers  are  %nrltten  In  specific  %«ys  to  define  their  type : 

1  .  Liteger:  2 

2  .  Real :  2  .  or  2  JO 

3  .  Standard  notation  2  .78  x  10^ :  2  .78  E-K)5 


(double  precision)  2  .78  DfOS 


I.  mtcDucrKii 


The  polygaonu  functions  appear  In  nuaerous  oatheaatlcal  expreaalons 
related  to  the  evaluation  of  streaoes  In  tapered  sections  of  gun  tubes  .  The 
Ballistic  Research  laboratory  (BRI^  has  need  of  special  function  subroutines 
such  as  the  one  described  In  this  report  for  use  In  large  nuaerlcal  analysis 
computer  codes  for  stress  simulations  and  experimental  data  evaluation .  There 
are  some  specific  requirements  placed  on  such  subroutines.  In  particular, 
memory  requirements  should  be  minimal,  since  part  of  the  code  may  be  run  on 
small  computers  In  the  pre-  or  post-processing  of  the  data  .  For  this  same 
reason,  the  subroutine  should  be  as  machine  Independent  as  possible.  That  Is, 
It  should  contain  no  "special  tricks"  related  to  a  specific  machine.  Finally, 
the  code  must  achieve  the  greatest  precision  possible,  since  the  output  Is 
merely  a  component  of  other  computations  . 

The  polygamma  function  subroutine  described  In  this  report  satisfies 
these  conditions  over  a  broad  range  of  application  requirements  .  It  Is 
designed  for  complex  arguments  In  double-precision;  the  subroutine  Is  written 
In  FCBITRAN  IV  and  does  not  require  Implicit  complex  arithmetic  .  Examples  run 
on  the  CDC  7600  computer  at  BRL  have  provided  a  precision  of  20-22  significant 
digits  .  Examples  run  on  an  HP-1000  minicomputer  at  BRL  have  provided  a 
precision  of  13-15  significant  digits. 

II.  nPOT  iMD  OOTPDT  VARIABUS 

The  subroutine  statement  Is 

SUBROUTINE  P0LYG<X,  Y,  N,  P0LYR,P0LYI,IERR)  . 

The  Input  variables  are  X,  Y,  and  N.  X  and  Y  are  the  real  and  imaginary 
parts,  respectively,  of  the  complex  argument  >X-flY  at  which  the  function  Is 
to  be  evaluated  .  X  and  Y  are  both  double-precision  real  variables  .  N  Is  an 
Integer  variable  which  specifies  which  polygamma  function  Is  to  be 
evaluated  .  Specifically, 

N  -  0  specifies  computation  of  the  Psi  function; 

N  -  1 ,2 ,3 ,..  .specifies  computation  of  the  derivative  of  the  Psl  function 
of  order  1,  2,  3,...,  respectively. 

NOTE :  If  >X+1Y  Is  In  the  left  half  plane ,  then  N  must  be  less  than  5  . 

The  output  variables  are  P0LYR,P0LYI,  and  lERR .  POLYR  and  POLYI  are 
double-precision  real  variables  which  are  the  real  and  Imaginary  parts, 
respectively,  of  the  requested  function,  lERR  Is  an  Integer  variable 
specifying  errors  detected  by  the  subroutine  .  Specifically,  If 

lERR  -  0,  no  errors  occurred; 

>  1,  negative  N  was  requested; 

■  2,  Input  argument  X-flY  tms  zero  or  a  negative  Integer; 

■  3,  N>4  and  X-flY  was  In  the  left  half  plane. 
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Appropriate  error  Beaeagee  are  printed  by  the  aubroutlne  to  aeeeBpany  the 
nonzero  valuea  of  lERR* 


III.  MRioD  or  oaaavavm 

A.  OeflnltloM.  The  gaaaa  function.  FCa).  le  defined  by* 


r(Z) 


llB 

n-** 


nin 


Z 


Z(Z+1)  ...  (Z4n) 


Z  *  0. 


(1) 


The  Pal  functlea.  ^(Z)i  la  defined  by 


it»(z)  -  r'(z)/r(z). 


(2) 


4i(z)  la  the  polygaaM  function  of  order  zero, 
polygaama  function  ''(Z)  la  defined  by 


In  general .  the  nth-order 


^(n)(Z)  •  (4.(2))  ,  (3) 

if 


where  n  la  a  poaltlve  Integer. 

The  method  of  coaputatlon  used 
asymptotic  series  for  4*^"  vZ).  If 
formulas  arc  used. 

B.  Asymptotic  Series.  For  values  of  Z  not  on  the  negative  real  axis, 


by  the  subroutine  la  evaluation  of  the 
the  magnitude  of  Z  la  small,  recursion 


nl 


2Z 


,n+l 


I 

k-1 


(2k-bn-l)! 

(2k)IZ^’‘'^ 


(4) 


*Ail  of  the  fotmiTae  used  in  thie  report  may  be  found  in  the  Handbook  of 
Mathematiaal  Funotione  of  the  national  Bureau  of  Standards. 


as  Z-»*  ,  where  the  symbol  I  denotes  the  factorial  operator,  and  the  821^  are 
the  Bernoulli  numbers,  defined  later  In  this  section.  In  the  case  n  ■  0, 


HZ) 


0* 

I 

n-1 


2n 


2nZ 


2n  * 


(5) 


C.  Recurrence  and  Reflection  Formulas.  Through  trial  computer  runs  It  was 
found  that  the  asymptotic  series  was  most  effective  for  values  of  Z  such  that 
jZ  |>10.  Vftten  |Z  |<10  the  subroutine  uses  the  recurrence  formula 


♦^"\z+l)  -  *^"^(2)  +  (-l)"n!  Z"““^ 


(6) 


of 


jl^en  Re  {z}  is  negative,  and  Z  Is  n 
(j)'“^(Z)  are  defined  In  terms  of 


ot  an  Integer,  the  values 

(1-Z)  by  use  of  the  reflection  formula 


♦^"\z)-(-l)%^"^(l-Z)  -  u  ^  cot  (nZ). 

dZ 


(7) 


It  was  decided  that  the  analytic  expression  for  the  first  four  derivatives  of 
the  cotangent  function  would  be  programmed  directly,  hence  the  restriction  on 
N. 

D.  Programming  Methods.  The  line  numbers  given  In  this  section  refer  to  the 
program  listing  In  Appendix  A.  For  operation  on  various  computers,  the  value 
of  PI  In  line  11  can  be  increased  or  decreased  In  precision.  It  should  be  as 
precise  as  a  particular  computer  will  allow.  Lines  12-30  define  the  Bernoulli 
numbers  B2,  B^,...B2g  for  use  In  the  asymptotic  series.  These  numbers  are 
expressed  as  quotients,  with  exact  numerators  and  denominators.  This  allows 
them  to  be  computed  to  the  full  precision  of  the  computer.  Clearly,  then,  no 
more  than  38  terms  of  the  asymptotic  series  are  used.  In  fact.  In  all  ranges 
of  values  of  the  argtiment  used  for  testing,  26  terms  was  the  maximum  ever 
required  to  achieve  optimal  precision.  Line  31  sets  this  maximum  number  of 
terms  to  38. 
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Line  32  sets  the  cutoff  value  below  which  recurrence  oust  be  used.  Lines 
33-50  check  for  Input  errors.  If  |Z|  Is  within  l.D-30  of  sero  or  a  negative 
Integer,  an  error  code  IERR-2  Is  set.  This  value  (line  33)  nay  be  changed  to 
suit  a  particular  computer's  range.  Lines  41-46  convert  Z  from  the  left  half 
plane  to  the  right  half  plane,  If  necessary.  Lines  47-50  convert  Z  to  the 
first  quadrant.  If  necessary,  to  simplify  the  computations.  The  sign  of  the 
imaginary  part  of  Z  is  saved  In  the  variable  SGNI  for  use  In  line  58. 

Evaluation  of  the  appropriate  function  Is  accomplished  In  lines  51-67. 

If  |Z|  Is  less  than  the  cutoff  value  set  In  line  32,  then  SUBROUTINE  RECUR  Is 
called  to  Implement  the  recurrence  Formula  (6).  If  N>0,  SUBROUTINE  PGAM  Is 
called  to  evaluate  the  polygamma  function.  If  N>0,  SUBROUTINE  PSI  is  called 
to  evaluate  the  Psl  function.  In  lines  57  and  58,  the  real  (ADDR)  and 
imaginary  (AOUI)  parts.  If  any,  due  to  recurrence  are  added  to  the  function 
values  POLYR  and  POLYI,  respectively.  In  line  58,  the  sign  of  the  imaginary 
part  Is  corrected  for  the  quadrant  In  which  Z  lies.  Lines  60-67  handle  the 
reflection  from  the  left  half  plane,  using  SUBROUTINE  COTAN  to  evaluate  the 
appropriate  derivative  of  the  cotangent  function. 

Each  of  the  subroutines  will  be  discussed  In  detail.  SUBROUTINE  RECUR 
(lines  86-111)  Implements  Formula  (6).  For  N>0,  lines  94-97  compute  N! . 

Lines  98-106  sum  the  powers  of  Z  computed  In  SUBROUTINE  CPOWR  (lines  193- 
206).  In  lines  108  and  109,  the  sums  of  the  real  and  Imaginary  parts  of  the 
powers  of  Z,  ADDR  and  ADDI,  respectively,  are  multiplied  by  NI  and  SGN,  which 
la  (-1)". 

SUBROUTINE  PGAM  (lines  112-165)  Implements  Formula  (4)  N>0.  Lines  119- 
123  compute  NI  and  (N-l)l.  Lines  124-130  compute  Z“^  and  V2  multiply 

them  by  (N-1)!  and  N! ,  respectively,  and  sum  these  terms.  The  remainder  of 
the  subroutine  sums  the  terms  of  the  series,  ending  with  NTERMS  or  when  the 
magnitude  of  the  Individual  terms  stops  decreasing. 

SUBROUTINE  PSI  (lines  166-192)  implements  Formula  (5).  It  uses  the  same 
technique  for  summing  the  asymptotic  series  as  was  used  in  SUBROUTINE  PGAM. 
SUBROUTINE  COTAN  (lines  207-256)  evaluate  the  appropriate  derivative  of  the 
cotangent  function,  according  to  the  formulas 


dz 


cot(irz) 


-  Tr(  1  +  cot^(irz)) 


(8) 


2  3 

— z-  cot(az)  ■  2ir^  (cot(irz)  +  cot'*(irz))  , 
dz*^ 


(9) 


(10) 


T  cot(itz)  ■  -  2w^(l  +  4  cot^(itz)  +  3  cot^(ir*))  , 
dz^ 


cot(iiz)  “  8it^(2  cot(iiz)  +  5cot^(irz)  +  3cot^(*z))  ,  (11) 

dz^ 


where 


cot  (irz) 


48ln(zx)co8(zx) 
e^^'^+e  ^’'^+48ln^xx-2 


+ 


ie 


-2xy_  2xy 


I 


2nyj^“2xy,.  .2  - 

i  +e  '+48ln  irx-2 


(12) 


for  z  «  x+ly. 


IV.  OMCUISKHIS 

The  eubroutine  deecrlbed  In  thl8  report  providee  valuee  of  the  polyganma 
functlona  for  complex  argument.  Although  precialon  will  vary  from  one 
computer  to  another,  the  eubroutine  hae  achieved  22-dlglt  precision  on  the  CDC 
7600. 
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SUBROUTINE  POLVG<X,Y,N,POLYR,POLYI . lERR) 

C 

C  ASYMPTOTIC  SERIES  SUBROUTINE  FOR 
C  THE  POtYGAMMA  FUNCTION 
C 

C  J.MALBERT,  OCT,  1981 

C 

IMPLICIT  DOUBLE  PRECISION  CA-H.O-Z) 

DIMENSION  B(52) 

PI -3. 1415926535897932384D0 

eC2)«1.00/6.00 

6C4)--1.DO/30.00 

BC6)*1. 00/42. 00 

B<8)--1 .00/30.00 

BC10)-5. 00/66.00 

8 (12) >-691. 00/2730. DO 

B(14)>7.DO/6.00 

B( 16) >-361 7. 00/510. DO 

B( 18) -43867. 00/ 798. DO 

8 (20)- -17461 1.00/330. DO 

8(22) -854513. 00/138. DO 

B (24) --236364 091. 00/2730. DO 

8(26) -8553103. DO/6. DO 

B (28)- -23749461 029. 00/870. DO 

B(30)-8615841276005. DO/14322. 00 

B(32) --7709321041217. 00/510. DO 

B(34) -2577687858367. 00/6. DO 

B (36) --26315271553053477373.00/ 1919190. DO 

B(38)-2929993913841559. 00/6.00 

NTERMS-38 

CUTOFF-10.00 

ZERO-1. 0-30 

ZR-X 

ZI  -V 

IF(N.LT.O)  GOTO  100 

IF(DABS(Zl).GT.ZERO.OR.ZR.GT.ZERO)  GOTO  1 
INT-ZP 

TEST-OBLE<FLOAT(INT)) 

IF(OABS(TEST-ZR).LE.ZERO)  GOTO  200 
1  REFLEC--1.D0 

IF(ZR.GE.O.DO)  GOTO  5 

IF(N.GT.4)  GOTO  300 

REFLEC-l.DO 

ZR-l.OO-ZR 

ZI--ZI 

5  SGNI-1.00 

IF(2I .GE. 0.00)  GOTO  10 

SGNl-l.OO 

ZI-DABS<ZI) 

10  AOOR-O.OO 

ADDI-O.DO 

W-OSQRT(ZR»ZR-»ZI*ZI  ) 

IF(W. LT. CUTOFF)  CALL  RECUR (W , ZR . Z I , N , CUTOFF , AOOR , ADO  I ) 

15 


fflKmiMa  Fiat  suMi-tior  num 


0055 

IFfN.'iT.O)  CALL  PGAMCZR ,ZI ,N ,NTERMS , B .POLYR ,POLV I) 

0056 

IF(N.EQ.O)  CALL  PS  I (ZR . Z I , HTERMS , B . POLYR ,POL Y I ) 

0057 

POLYR-POLYR+AOOR 

0058 

POLYI-SGHI*<POL VI *80017 

0059 

IFCREFLEC.lt. 0.00)  60T0  999 

0060 

IFCMOD(H,2) .EQ.O)  GOTO  15 

0061 

P0LVR-*P0LYR 

0062 

POLY  I --POLY  I 

0063 

15 

ZR-X 

0064 

21-Y 

0065 

CALL  COTANCZR.ZI .N.P! , COTDR ,COTD I ) 

0066 

POLYR-POLYR-COTOP 

0067 

POLYl-POLYI-COTOI 

0068 

GOTO  999 

0069 

100 

IERR-1 

0070 

ImRITECI  ,1000) 

0071 

GOTO  999 

0072 

200 

lEPR-2 

0073 

WRITECl ,2000) 

0074 

GOTO  999 

0075 

300 

IERR-3 

0076 

WRITECl ,3000) 

0077 

999 

RETURN 

0078 

1000 

FORMATC •**ERROR  FROM  POLVGAMMA  SUBROUT  I NE««* ,/ , 

0079 

*  *  NEGATIVE  N  REQUESTED*) 

0080 

2000 

FOPMATC  *«ERROR  FROM  POLYGAMMA  SUBROUT  I  NEK** ,/ , 

0081 

«  •  INPUT  ARGUMENT  MAS  ZERO  OR  A  NEGATIVE  INTEGER*) 

0082 

3000 

FORMAT <•  »»ERROR  FROM  POLYGAMMA  SUBROUT  I NE»»* , 

0083 

*  "INPUT  ARGUMENT  WAS  IN  THE  LEFT  HALF  PLANE,*,/, 

0084 

»  *  WITH  REQUESTED  ORDER  N  GREATER  THAN  4*) 

0085 

END 

0086 

SUBROUTINE  RECUR  CM, ZR ,21 , N , CUTOFF , ADOR , ADD i ) 

0087 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z) 

0038 

NTIMES-CUTOFF-W 

0089 

IFCNTIMES.lt. 1)  NTIMES-1 

0090 

SGN=-1.00 

0091 

FACT-1 .00 

0092 

IFCN.LE.O)  GOTO  10 

0093 

IF(MOD<H,2) .HE.O)  SGH-1.00 

0094 

EN-DBLE(FLOAT(N)) 

0095 

DO  5  I-1,N 

0096 

FACT«FACT*EN 

0097 

5 

EN-EN-1.00 

0098 

10 

ZRADD-O.DO 

0  099 

tC-N*l 

0100 

00  15  l-l,NTIMES 

0101 

CALL  CPOWR(ZR*ZRAOO.ZI ,ZPR,ZPI ,K) 

0102 

U-ZPR*ZPR*2PI«ZPI 

0103 

AODR-ADOR*ZPR/U 

0104 

ADOI-AODI-ZPI/U 

0105 

ZRADD-DBLECFLOATCD) 

0106 

15 

CONTINUE 

0107 

ZR-ZR*2RAD0 

0108 

ADC>R-ADDR»SGN«FACT 

16 


0109 

ADD I -ADD  I vSCHaF ACT 

0110 

RETURN 

0111 

END 

0112 

SUBROUTINE  PGAH<ZR,Z1 .N.NTERNS.B.POLVR.POLYI) 

0113 

IMPLICIT  DOUBLE  PREC’SION  CA-H,0-Z) 

0114 

DIMENSION  B<S2> 

0115 

PROD- 1. DO 

one 

HSTOP-N-1 

0117 

EN-l.DO 

one 

IFtNSTOP.LT.n  GOTO  25 

one 

DO  20  I-l,NSTOP 

0120 

PRQD-PROOBEN 

0121 

20 

EN-EN-l.DO 

0122 

25 

ENMl-PROD 

0123 

PROD-PROO»EN 

0124 

CALL  CP0MR(ZR,2I.ZNR,ZNI,N} 

0125 

U-ZNRaZNR«2NI«ZNI 

0126 

ZNRP 1 - ZNRaZR - ZN I »Z I 

0127 

ZN I P 1 -ZNRaZ I «ZN I «ZR 

012S 

U-2 .  DOXZNRPIVZNRPI  «ZN  I  PiaZN  1  PI ) 

0129 

P0LYR-PR0D«ZNRP1/U-ENM1«ZNR/M 

0130 

POLVl --PR00«ZNIP1/U-ENM1«ZNI /M 

0131 

EK-l.OO 

0132 

SUMR-O.OO 

0133 

SUMI-0.00 

0134 

FACT-ENMl 

0135 

BASE-EN 

0136 

DFACT-l.DO 

0137 

ZRLST- 1.030 

0138 

ZILST-1.030 

0139 

DO  100  K-2,NTeRMS,2 

0140 

FACT-FACT»BASE 

0141 

BASE-BASE*!. 00 

0142 

FACT-FACT»BASE 

0143 

BASE-BASE*!. DO 

0144 

OFACT-OFACT»EK 

0145 

EK-EK*1.00 

0146 

OFACT-OFACT«EK 

0147 

EK-EK*1.00 

0148 

COEF-BCK)«F ACT/OF ACT 

0149 

CALL  CPOMR(ZR>ZI.ZNR,ZNI,K*N> 

0150 

M-ZNRaZNR*ZNlaZNI 

0151 

ZNR-COEF«ZNR/U 

0152 

ZNI--COeF«ZHI/M 

0153 

IF(OABS(ZNR>.GT. ZRLST. OR. DABSCZND.GT.ZILST)  GOTO  110 

0154 

ZRLST-DABSCZNR> 

0155 

ZILST-DABSCZNI) 

0156 

SUMR-SUMR*ZNR 

0157 

SUMI-SUMI*ZNI 

0158 

100 

CONTINUE 

0159 

no 

POLYR-POLYR*SOMR 

0160 

POLYI-POLYI*SUMI 

0161 

IF(MOD(N.2>.NE.O}  GOTO  200 

0162 

POLYR--POLVR 
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0163 

POLYI--POLVI 

0164 

200 

RETURN 

0165 

END 

0166 

SUBROUTINE  PSIlZR.ZI ,NTERMS.B,POLVR.POLVI) 

0167 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

0168 

DIMENSION  B<52) 

0169 

M-ZR»ZR*ZI«ZI 

0170 

irCH.LT. 1.0-30)  GOTO  99 

0171 

POLVR-DLOG(DSQRT<M)>-.5DO»ZR/M 

0172 

POLY I -OATAN  <Z I /ZR) « . 5D0aZ I /H 

0173 

SUMR-O.OO 

0174 

SUMI-O.DO 

0175 

ZRLST- 1.030 

0176 

Z1LST-1.030 

0177 

DO  10  K-2.NTERMS,2 

0178 

COEr>B(K)/OBLECFLOAT(K)) 

0179 

CALL  CPOMR(ZR,ZI,ZNR,ZNI,K) 

0180 

M-ZNRVZNR«ZNI«ZNI 

0181 

ZNR-COEF«ZNR/U 

0182 

ZNI*-COEF»ZNI/U 

0183 

I F  COABS(ZNR) .GT . ZRLST. OR. OABS(ZNI ) . GT. ZI LST) 

0184 

ZRLST>OABS(ZNR) 

0185 

ZILST-OABS(ZNi) 

0186 

SUMR-SUMR«ZNR 

0187 

SUM1-SUMI«ZNI 

0188 

10 

CONTINUE 

0189 

20 

POLYR-POLYR-SUMR 

0190 

POLYI-POLYI-SUMI 

0191 

99 

RETURN 

0192 

END 

0193 

SUBROUTINE  CPOWR<X,Y,U,V,NTIMES) 

0194 

DOUBLE  PRECISION  X,Y,U,V,H,S,T 

0195 

S-X 

0196 

T-Y 

0137 

U-S 

0198 

V-T 

0199 

IFCNTIMES.LE. 1)  GOTO  20 

0200 

NTMl-NTlMES-1 

0201 

DO  10  I-1,NTM1 

0202 

W-S*U-T*V 

0203 

V-T«U«S»V 

0204 

10 

U>H 

0205 

20 

RETURN 

0206 

END 

0207 

SUBROUT  I NE  COTANIZR . Z I , N , P I , COTDR , COTD I ) 

0208 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z} 

0209 

HI-PI<iZR 

0210 

HCOS>OCOS(MI) 

0211 

HSIN>DSIN(HI) 

0212 

IFCZI.EQ.0.00)  GOTO  1 

0213 

UINV-DEXP(P1«Z1) 

0214 

U-l.OO/UINV 

0215 

DENOM-UvU«UINVaUINV«4.DO«HSIN«MSIN-2.00 

0216 

CTPZR-4 .  DOMIS 1  N«HCOS/DENOM 
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0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 


CTPZ I - (U>U-U I NVMU I NV> / DENON 
GOTO  5 

1  CTPZR-UCOS/MSIN 

CTPZI-0.00 

5  IF(N.EQ.0>  GOTO  60 

IFCN.GT.4>  GOTO  50 
GOTO  <10, 20, 30, 40), N 

10  CALL  CPQI4RCCTPZR,CTPZI,COTDR,COTOI ,2) 
COTOR--PI«(1.DO«COTDR) 

COTOl— PIBCOTOI 
GOTO  70 

20  CALL  CPOHR<CTPZR,CTPZI ,COTDR,COTDI ,3) 

C0EF-2.00«P1«PI 
COTDR-COEF«<CTPZR*COTDR> 

COTD I -COEF«<CTPZ I ♦COTD I > 

GOTO  70 

30  CALL  CPOMR<CTPZR,CTPZI,TEMPR,TEHPI,2) 

CALL  CPOMR<TENPR.TENPI ,COTDR,COTDI ,2) 
COEF--2.00aPI«PlaPI 

COTDR>COEF»<l.DO«4.DOaTEHPR«3.DO»COTDR) 

COTS  I -C0EFV<4. 00«TEHP t  «3. DOaCOTO I > 

GOTO  70 

40  CALL  CPOMR<CTPZR,CTPZI ,TEHPR,TEMPI ,3) 

CALL  CPOMR<CTPZR,CTPZI,COTDR,COTDI ,2> 
U-COTDR«TEI4PR-COTD  I  aTEHP  I 
COTR I -COTD I aTEMPR^COTDRaTEMP I 
COTDR-U 

COEF -8 . DO»P I »P J »P I *P 1 

COTOR-COEFa(2.DOaCTPZR*5.DO»TEnPR*3.00aCOTDR> 
COTOl -C0EFa(2. OOaCTPZI «5. DOaTEHPl «3. OOaCOTOI > 
GOTO  70 

50  COTOR-O.OO 

COTD I -0.00 
GOTO  70 

60  COTOR-CTPZR 

COTOI-CTPZI 
70  COTDR-PlaCOTOR 

COTOI-PlaCOTOI 
RETURN 
END 
ENDS 
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